# ------------------------------------------------------------------------------
# Produces reg tables for symptom salience, representativeness
# Author: Cassidy Shubatt <cshubatt@gmail.com>
# To run: bash 02_symptom_salience.sh
# ------------------------------------------------------------------------------

# Libraries --------------------------------------------------------------------
message("Loading libraries...")
library(here)
library(yaml)
library(tidyverse)
library(glue)
library(testit) # assert()
library(broom) # tidy()
library(data.table) # setnames()
library(lfe) # felm
library(xtable)
library(stargazer)

temp <- here("code", "07_other_physician_errors", "temp")
u <- modules::use(here("lib", "util.R"))

# Load Data --------------------------------------------------------------------
message("Loading data...")
overnight_lab <- ""
paths <- read_yaml(here::here("lib", "filepaths.yml"))
cc <- readRDS(paths$analysis$cc) %>% select(-cc_arrest_cardiac)
cc_labels <- read_csv(here("lib", "cc_labels.csv"))
cc_freq <- readRDS(file.path(temp, "cc_frequency_df.rds"))

cohort <- readRDS(glue(paths$analysis$test_cohort)) %>%
  u$safe_left_join(cc) %>%
  filter(!exclude)

# Prep  ------------------------------------------------------------------------
message("Preparing data for regressions...")
model_df <- cohort %>%
  mutate(yhat_full = p__ensemble__stent_or_cabg_010_day__tested)

# Individual Symptom Regressions -----------------------------------------------
message("Regressing Y, T on yhat and individual symptoms...")
cc_topx <- cc_freq$chief_complaint[1:10]
outcomes <- c("test_010_day", "stent_or_cabg_010_day")

# reg tagble row labs
cc_labels <- cc_labels %>%
  filter(chief_complaint %in% cc_topx)
covars <- c("yhat_full", cc_labels$chief_complaint)
covar_labs <- c("Predicted Risk", cc_labels$cc_label)
names(covar_labs) <- covars
message("Covariate Labeling Scheme: ")
print(covar_labs)

reg_list <- list()
for(outcome in outcomes){
  form <- reformulate(
    response = outcome,
    termlabels = c("yhat_full", cc_topx, "1 | 0 | 0 | ptid")
  )
  reg_df <- copy(model_df)
  if(outcome == "stent_or_cabg_010_day"){
    reg_df <- filter(reg_df, test_010_day)
  }else if(outcome == "macetrop_030_pos"){
    reg_df <- filter(reg_df, !test_010_day)
  }

  fit <- felm(form, reg_df)
  reg_list <- c(reg_list, list(fit))
}

reg_tbl <- stargazer(
  reg_list, omit.stat = c("f", "ser", "ll", "aic", "adj.rsq"),
  dep.var.caption = "", order = names(covar_labs),
  covariate.labels = covar_labs
)
save_fp <- file.path(temp, "symptom_salience_regressions.tex")
write(reg_tbl, save_fp, append = FALSE)

message("Done.")
